### REPLICATION CODE FOR DUNHAM, LIEBERMAN, AND SNELL (2016)
### Analyses and figures in Appendix S3

setwd("~/Dropbox/DLS_racerisk_PLOSms/replication_data") ##Set your own directory

# required R packages
library ("foreign")
library("ggplot2")
library("plyr")
library("Rmisc")
library("reshape2")
library("effsize")

# load data
risk <- read.dta("DLS_riskdata_tot.dta") 


## Appendix S3 Section 1 From Here

# Correlations between Time 1 and Time 2 identification scores
cor(risk$id_AA,risk$id_AA2,use='complete.obs')
cor(risk$id_white,risk$id_white2,use='complete.obs')
cor(risk$id_dem,risk$id_dem2,use='complete.obs')
cor(risk$id_old,risk$id_old2,use='complete.obs')
cor(risk$id_poor,risk$id_poor2,use='complete.obs')
cor(risk$id_rich,risk$id_rich2,use='complete.obs')
cor(risk$id_rep,risk$id_rep2,use='complete.obs')
cor(risk$id_young,risk$id_young2,use='complete.obs')

# statistical significance of lowest correlation from above
cor.test(risk$id_young,risk$id_young2)

# average correlation of above correlations
(cor(risk$id_AA,risk$id_AA2,use='complete.obs')+cor(risk$id_white,risk$id_white2,use='complete.obs')+
  cor(risk$id_dem,risk$id_dem2,use='complete.obs')+cor(risk$id_old,risk$id_old2,use='complete.obs')+
  cor(risk$id_poor,risk$id_poor2,use='complete.obs')+cor(risk$id_rich,risk$id_rich2,use='complete.obs')+
  cor(risk$id_rep,risk$id_rep2,use='complete.obs')+cor(risk$id_young,risk$id_young2,use='complete.obs'))/8

# create summary statistics for the two main race identification items (Time 1)
id_AA <- summarySE(risk,measurevar='id_AA',groupvars='black')
names(id_AA) <- c('race','N','identity','sd','se','ci')
id_AA$type <- 'AfAm'
id_white <- summarySE(risk,measurevar='id_white',groupvars='black')
names(id_white) <- c('race','N','identity','sd','se','ci')
id_white$type <- 'White'
raceident <- rbind(id_AA,id_white)

# Figure A from Appendix S2
ggplot(raceident, aes(x=type, y=identity, fill=race)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=identity-ci, ymax=identity+ci),
                width=.2,                   
                position=position_dodge(.9)) + 
  scale_fill_hue(name="Participant Race",
                 labels=c("White", "Af-Am")) +
  xlab("Target of Identification") + ylab("Degree of Identification (0 to 100)") +
  ggtitle("Pre-experiment Racial Identification, \nby participant race and target of \nidentification (+ 95% CIs)") +
  expand_limits(y=c(0,100)) +
  scale_y_continuous(breaks=0:100*25) +
  scale_x_discrete(labels=c("Af-Am", "White"))

# compute means and SDs for identification items, by respondent race
identity <- risk[c(57,74:81)]
idents <- melt(identity,id.vars=c("black"))
id_means <- summarySE(idents,measurevar='value',groupvars=c('black','variable')) # this dataframe has means and SDs from Table 1 (Appendix S3)

# compute effect sizes for race difference, as reported in Table 1 (Appendix S3)
cohen.d(risk$id_old,risk$black)
cohen.d(risk$id_AA,risk$black)
cohen.d(risk$id_white,risk$black)
cohen.d(risk$id_dem,risk$black)
cohen.d(risk$id_rep,risk$black)
cohen.d(risk$id_young,risk$black)
cohen.d(risk$id_rich,risk$black)
cohen.d(risk$id_poor,risk$black)

# p-value associated with smallest differences
t.test(risk$id_old~risk$black) 
t.test(risk$id_rich~risk$black) 

## End Appendix S3 Section 1

## Appendix S3 Section 2 From Here

# rename critical vars for my analysis with meaningful names
risk=rename(risk, c("bw1"="sepWorlds", "bw2"="diffNeeds", "bw3"="distComms","bw4"="decrBounds"))
# convert to numeric
risk$sepWorlds <- as.numeric(risk$sepWorlds)
risk$diffNeeds <- as.numeric(risk$diffNeeds)
risk$distComms <- as.numeric(risk$distComms)
risk$decrBounds <- as.numeric(risk$decrBounds)

# correlation between racial polarization items
cor(risk[85:88])

# compute means and SDs for polarization items, by respondent race
bw <- risk[c(57,85:88)]
bws <- melt(bw,id.vars="black")
bw_means <- summarySE(bws,measurevar='value',groupvars=c('black','variable'))

# Produce Figure B from Appendix S3
ggplot(bw_means, aes(x=variable, y=value, fill=black)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=value-ci, ymax=value+ci),
                width=.2,                   
                position=position_dodge(.9)) + 
  scale_fill_hue(name="Participant Race",
                 labels=c("White", "Af-Am")) +
  xlab("Item") + ylab("Degree of Endorsement (1 to 5)") +
  ggtitle("Means and 95% CIs for Racial Polarization Items \nby Respondent Race") +
  expand_limits(y=c(1,5)) +
  scale_y_continuous(breaks=1:5*1) +
  scale_x_discrete(labels=c("Separate\nworlds", "Different\nneeds", "Distinct\ncommunities", "Decrease\nboundaries"))

# tests of race differences
t.test(risk$sepWorlds~risk$black) 
t.test(risk$diffNeeds~risk$black) 
t.test(risk$distComms~risk$black) 
t.test(risk$decrBounds~risk$black)

# effect sizes associated with race differences
cohen.d(risk$sepWorlds,risk$black)
cohen.d(risk$diffNeeds,risk$black)
cohen.d(risk$distComms,risk$black)
cohen.d(risk$decrBounds,risk$black)

## End Appendix S3 Section 2




